import numpy as np


def create(n):
    """
    Create a n-dim Extended Rosenbrock function

    :param n: dimensionality
    :return: x0, f, g, G
    """
    if n % 2 != 0:
        raise ValueError("n must be even")

    x0 = np.zeros(n, dtype=np.float64)
    x0[0::2] = -1.2
    x0[1::2] = 1.0

    def f(x: np.ndarray):
        assert x.ndim == 1
        assert x.shape[0] == n
        x1 = x[0::2]
        x2 = x[1::2]
        f1 = 100 * (x2 - x1 ** 2) ** 2
        f2 = (1 - x1) ** 2
        return np.sum(f1 + f2)

    def g(x: np.ndarray):
        assert x.ndim == 1
        assert x.shape[0] == n
        x1 = x[0::2]
        x2 = x[1::2]
        grad = np.zeros(n, dtype=np.float64)
        grad[0::2] = -400 * (x2 - x1 ** 2) * x1 - 2 * (1 - x1)
        grad[1::2] = 200 * (x2 - x1 ** 2)
        return grad

    def h(x: np.ndarray):
        assert x.ndim == 1
        assert x.shape[0] == n
        x1 = x[0::2]
        x2 = x[1::2]
        h11 = -400 * (x2 - 3 * x1 ** 2) + 2
        h12 = -400 * x1
        base = 2 * np.arange(0, n // 2, dtype=np.int32)
        he = np.zeros((n, n), dtype=np.float64)
        he[base + 0, base + 0] = h11
        he[base + 1, base + 0] = h12
        he[base + 0, base + 1] = h12
        he[base + 1, base + 1] = 200
        return he

    return x0, f, g, h
